%Edited 11/15/20 to just output the numerical quantification data, (normalized)
%% This script quantifies the protein distribution inside the Lymph node
% depending on the distance to the edge of the Lymph node.
% 08/10/2020
clear all 
close all
mkdir Graphs
%% Variables
px_size=0.722; %pixel size
step = 10; %the step that we use is the ring radius to measure intensity vs distance to edge
step_sq = 70; %the side of the square to measure the local distribution of intensity
%% Load image:
for i=1:3
image(:,:,i)=imread('a',i);
end
imshow(image);
close;

%% draw region of interest
BW=roipoly(image);
%imshow(BW);
%close;

%% Intensity vs distance from the edge
% We define the radius of the ring where the average intensity will be
% measured

[Dista,Linear] = bwdist(~BW);
Dista = round(Dista); % Distances from the edge.

% just to see that the distance is properly done we plot Distance
imagesc(Dista);
colorbar;
%print(['distance','.png'],'-dpng','-r150');
close;

%We quantify the intesity per ring

%We quantify it per channel
for z=1:1:length(image(1,1,:))
mean_image=zeros(size(Dista));
mean_per_region=zeros(size(image));        
j=1;   
for k = 1:step:max(Dista(:))
    dist_reg(j,1)=k;
    d_temp = Dista(:,:) < k+step & Dista(:,:) > k-1;
    int_ring= image(:,:,z).*uint8(d_temp);
    mean_int_ring(j)=mean(int_ring(find(int_ring)));
    mean_image=mean_image + mean_int_ring(j)*double(d_temp);
    j=j+1;
end
real_dista=dist_reg.*px_size;
mean_per_region(:,:,z)=mean_image(:,:);
mean_int_ring_channel(:,z)=mean_int_ring(:);
imagesc(mean_image); %plot mean value per region.
%{
colorbar;
title(['Mean intensity per region channel',num2str(z)])
print(['Graphs\mean int per region channel ',num2str(z),'.png'],'-dpng','-r150');
close;
plot(mean_int_ring); 
title(['Mean intensity vs distance to the lymph node edge',num2str(z)]);
xlabel('Distance to the edge');
ylabel('Mean intensity per region');
%print(['Graphs\graph mean int per region channel ',num2str(z),'.png'],'-dpng','-r150');
%}
    
end 

%plot (real_dista, mean_int_ring_channel(:,1),'mo','LineWidth',2);
%hold on;
%plot (real_dista, mean_int_ring_channel(:,2),'ro','LineWidth',2);
%hold on;
%plot (real_dista, mean_int_ring_channel(:,3),'ko','LineWidth',2);
%xlabel('Distance to the edge');
%ylabel('Mean intensity per region');
%legend('channel 1', 'channel 2', 'channel 3'); 
%print(['Graphs\graph mean int vs distance to edge.png'],'-dpng','-r150');

%%Normalize the values to their maximum values
for w=1:length(image(1,1,:))
    norm_int_per_region(:,w)= mean_int_ring_channel(:,w)./max(mean_int_ring_channel(:,w));
end    
%plot (real_dista, norm_int_per_region(:,1),'mo','LineWidth',2);
%hold on;
%plot (real_dista, norm_int_per_region(:,2),'ro','LineWidth',2);
%hold on;
%plot (real_dista, norm_int_per_region(:,3),'ko','LineWidth',2);
%xlabel('Distance to the edge');
%ylabel('Mean intensity per region');
%legend('channel 1', 'channel 2', 'channel 3'); 

%print(['Graphs\graph norm int vs distance to edge','.png'],'-dpng','-r150');
%close;

properties = [dist_reg real_dista norm_int_per_region mean_int_ring_channel];

%Displaying distance and normalized data
disp("distance")
disp(real_dista)
disp("Channel 1, raw")
disp(mean_int_ring_channel(:,1))
disp("Channel 2, raw")
disp(mean_int_ring_channel(:,2))
disp("Channel 3, raw")
disp(mean_int_ring_channel(:,3))
%save(['Properties.txt'],'properties','-ASCII');

%% Measure different channels local correlation

im_ch1 = zeros(numel(1:step_sq:size(image,1)),numel(1:step_sq:size(image,2)));
im_ch2 = zeros(numel(1:step_sq:size(image,1)),numel(1:step_sq:size(image,2)));
im_ch3 = zeros(numel(1:step_sq:size(image,1)),numel(1:step_sq:size(image,2)));

for g = 1:step_sq:size(image,1)-step_sq
    for n = 1:step_sq:size(image,2)-step_sq

    temp_ch1 = image(g:(g+step_sq-1),n:(n+step_sq-1),1);
    temp_ch2 = image(g:(g+step_sq-1),n:(n+step_sq-1),2);
    temp_ch3 = image(g:(g+step_sq-1),n:(n+step_sq-1),3);
    
    im_ch1((g-1)/step_sq+1,(n-1)/step_sq+1) = mean(temp_ch1(find(temp_ch1)));
    im_ch2((g-1)/step_sq+1,(n-1)/step_sq+1) = mean(temp_ch2(find(temp_ch2)));
    im_ch3((g-1)/step_sq+1,(n-1)/step_sq+1)=mean(temp_ch3(find(temp_ch3)));
    end
end

%{
subplot(2,3,1)
imagesc(im_ch1)
axis off
colormap(jet(256));
colorbar;
title('Mean local int in ch 1')
%caxis([0 15]);

subplot(2,3,2)
imagesc(im_ch2)
colormap(jet(256));
colorbar;
title('Mean local int in ch 2')
subplot(2,3,3)
imagesc(im_ch3)
colormap(jet(256));
colorbar;
title('Mean local int in ch 3')
%}

%% Normalize each channel to find local correlations between channels

%{
im_ch1_norm= im_ch1./max(im_ch1);
im_ch2_norm= im_ch2./max(im_ch2);
im_ch3_norm= im_ch3./max(im_ch3);

corr12=im_ch1_norm.*im_ch2_norm;
corr13=im_ch1_norm.*im_ch3_norm;
corr23=im_ch2_norm.*im_ch3_norm;

subplot(2,3,4)
imagesc(corr12)
axis off
colormap(jet(256));
colorbar;
title('correlation ch 1 and 2')
%caxis([0 15]);

subplot(2,3,5)
imagesc(corr13)
colormap(jet(256));
colorbar;
title('correlation ch 1 and 3')
subplot(2,3,6)
imagesc(corr23)
colormap(jet(256));
colorbar;
title('correlation ch 2 and 3')
print(['Graphs\intensity distribution','.png'],'-dpng','-r150');
close;
%}


    